function [fq, fmridata] = fmriquality(images, opts)
% fmriquality  - get some quality assurance for fMRI data
%
% FORMAT:       [fq, fmridata = ] fmriquality(images [, opts])
%
% Input fields:
%
%       images      list of images (or fMRI datatype, like FMR/VTC)
%       opts        optional settings
%        .motcor    perform motion-correction (and re-do stats, def: true)
%        .qasheet   flag, display quality assessment sheet (default: true)
%        .res       resolution (for motion detection, default: from file)
%        .tempffrq  temp filter frequency cut-off (as TRs, default: 80)
%        .tempfset  temp filter set, either of 'DC', {'Fourrier'}
%
% Output fields:
%
%       fq          complex struct, containing masks, time courses, etc.
%       fmridata    4-D data slab (motion corrected, if selected)

% Version:  v0.8a
% Build:    9110520
% Date:     Nov-05 2009, 8:46 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
   ((numel(images) ~= 1) || ...
    (~isBVQXfile(images, 'fmr') && ...
     ~isBVQXfile(images, 'vtc'))) && ...
   (~iscell(images) || ...
     numel(images) < 3)
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end
if nargin < 2 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'motcor') || ...
   ~islogical(opts.motcor) || ...
    numel(opts.motcor) ~= 1
    opts.motcor = true;
end
if ~isfield(opts, 'qasheet') || ...
   ~islogical(opts.qasheet) || ...
    numel(opts.qasheet) ~= 1
    opts.qasheet = true;
end
if ~isfield(opts, 'res') || ...
   ~isa(opts.res, 'double') || ...
    numel(opts.res) ~= 3 || ...
    any(isinf(opts.res) | isnan(opts.res) | opts.res <= 0 | opts.res > 16)
    opts.res = [];
else
    opts.res = opts.res(:)';
end
if ~isfield(opts, 'tempffrq') || ...
   ~isa(opts.tempffrq, 'double') || ...
    numel(opts.tempffrq) ~= 1 || ...
    isnan(opts.tempffrq) || ...
    opts.tempffrq < 8
    opts.tempffrq = 80;
end
if ~isfield(opts, 'tempfset') || ...
   ~ischar(opts.tempfset) || ...
    isempty(opts.tempfset) || ...
   ~any(lower(opts.tempfset(1)) == 'df')
    opts.tempfset = true;
else
    opts.tempfset = (lower(opts.tempfset(1)) == 'd');
end

% prepare output
fq = struct( ...
    'Dims', [], ...
    'Masks', struct, ...
    'Raw', struct, ...
    'TempFiltered', struct, ...
    'Quality', struct, ...
    'TC', struct, ...
    'MotCorr', struct);

% read in data (depending on type)
if numel(images) == 1
    if isBVQXfile(images, 'fmr')
        if images.FileVersion < 5 || ...
            images.DataStorageFormat < 2
            fmridata = images.Slice(1).STCData(:, :, :);
            fmridata(end, end, end, numel(images.Slice)) = 0;
            for sc = 2:size(fmridata, 4)
                fmridata(:, :, :, sc) = images.Slice(sc).STCData(:, :, :);
            end
            fmridata = double(permute(fmridata, [1, 2, 4, 3]));
        else
            fmridata = double(permute(images.Slice.STCData(:, :, :, :), [1, 2, 4, 3]));
        end
        if isempty(opts.res)
            opts.res = [ ...
                images.InplaneResolutionX, ...
                images.InplaneResolutionY, ...
                images.SliceThickness + images.GapThickness];
        end
    elseif isBVQXfile(images, 'vtc')
        fmridata = double(permute(images.VTCData(:, :, :, :), [4, 2, 3, 1]));
        if isempty(opts.res)
            opts.res = images.Resolution(1, [1, 1, 1]);
        end
    else
        % no other types yet
        error( ...
            'BVQXtools:NotYetSupported', ...
            'Currently only FMR and VTC BVQXfile types supported.' ...
        );
    end
elseif iscell(images)
    for sc = numel(images):-1:1
        if ischar(images{sc})
            try
                if ~isempty(regexpi(images{sc}(:)', '.*\.img$'))
                    images{sc} = [images{sc}(1:numel(images{sc})-3) 'hdr'];
                end
                images{sc} = BVQXfile(images{sc}(:)');
                if ~isBVQXfile(images{sc}, 'hdr')
                    error('NO_HDR');
                end
                ihdr = images{sc};
                images{sc} = ihdr.VoxelData(:, :, :, :);
                if isempty(opts.res)
                    opts.res = ihdr.ImgDim.PixSpacing(2:4);
                end
                ihdr.ClearObject;
            catch
                images(sc) = [];
            end
        elseif numel(images{sc}) == 1 && ...
            isBVQXfile(images{sc}, 'hdr')
            images{sc} = images{sc}.VoxelData(:, :, :, :);
            if isempty(opts.res)
                opts.res = ihdr.ImgDim.PixSpacing(2:4);
            end
        elseif ~isnumeric(images{sc}) || ...
            isempty(images{sc})
            images(sc) = [];
        end
    end
    try
        fmridata = double(cat(4, images{:}));
        [images{:}] = deal([]);
        images(1:end) = [];
    catch
        error( ...
            'BVQXtools:DimsMismatch', ...
            'Error concatenating images in 4th dimension.' ...
        );
    end
end
fmrisize = size(fmridata);
nslc = fmrisize(3);
nvol = fmrisize(4);
if nvol < 24
    error( ...
        'BVQXtools:TooFewImages', ...
        'Too few volumes in time series.' ...
    );
end
fq.Quality.Outliers.VolumeRatio = 0;
fq.Quality.Outliers.Volumes = zeros(nvol, 1);

% compute some basic images
fq.Dims = fmrisize;
fq.Raw.MeanImage = (1 / nvol) .* sum(fmridata, 4);
fq.Raw.StdevImage = std(fmridata, [], 4);

% get voxels that are likely to be background
mfmrimean = mean(fq.Raw.MeanImage(:));
mdfmrimean = median(fq.Raw.MeanImage(fq.Raw.MeanImage < mfmrimean));
fmribg = (fq.Raw.MeanImage < mdfmrimean);
fmrifg = (fq.Raw.MeanImage > mfmrimean);
fmrifg = fmrifg & dilate3d(erode3d(fmrifg));

% and remove those with more than 5% zero values
fmribg(sum(fmridata == 0, 4) > (0.05 * nvol)) = false;

% find biggest chunk
[cs, fmribg] = clustercoordsc(fmribg);
fmribg = (fmribg == maxpos(cs));
[cs, fmrifg] = clustercoordsc(fmrifg);
fmrifg = (fmrifg == maxpos(cs));
fq.Masks.Background = fmribg;
fq.Masks.Foreground = fmrifg;
nforeg = sum(fmrifg(:));

% create temporary mask to estimate smoothness over time
fmrisest = median(fq.Raw.MeanImage(fmrifg));
fmrisesd = 0.5 * std(fq.Raw.MeanImage(fmrifg));

% mask being within foreground where values between median +/- 0.5 * std
fmrisest = fmrifg & ...
    (fq.Raw.MeanImage > (fmrisest - fmrisesd)) & ...
    (fq.Raw.MeanImage < (fmrisest + fmrisesd));
fmrisest = smoothdata3(double(fmrisest), [2, 2, 2]) > 0.75;
[cs, fmrisest] = clustercoordsc(fmrisest);
fmrisest = (fmrisest == maxpos(cs));

% estimate smoothness (over time), iterate over vols
fq.TC.SmoothEst = zeros(nvol, 1);
for vc = 1:nvol
    
    % get volume data
    svol = fmridata(:, :, :, vc);
    
    % compute standard deviation within temp mask
    fmrisesd = std(svol(fmrisest));
    
    % actually smooth data
    svol = smoothdata3(svol, [2, 2, 2]);
    
    % and then recompute std and compare to unsmoothed value
    fq.TC.SmoothEst(vc) = std(svol(fmrisest)) ./ fmrisesd;
end
[trashvar, outlest] = windsorize(fq.TC.SmoothEst, 6, 3);
fq.Quality.Outliers.Volumes = fq.Quality.Outliers.Volumes + double(outlest(:));
fq.Quality.Outliers.VolumesSmoothEst = find(outlest);

% get some time courses
fq.TC.Global = lsqueeze(mean(mean(mean(fmridata))));
fq.TC.Slices = squeeze(mean(mean(fmridata)))';
fq.TC.Foreground = fq.TC.Global;
fq.TC.ForeSlices = fq.TC.Slices;
for vc = 1:nvol
    svol = fmridata(:, :, :, vc);
    fq.TC.Foreground(vc) = mean(svol(fmrifg));
    for sc = 1:nslc
        sslc = svol(:, :, sc);
        fq.TC.ForeSlices(vc, sc) = mean(sslc(fmrifg(:, :, sc)));
    end
end
fq.TC.ForeSlices(isnan(fq.TC.ForeSlices)) = 0;

% weigh slice timecourses by number of voxels
fq.TC.ForeSlicesWeighted = fq.TC.ForeSlices;
for sc = 1:nslc
    fq.TC.ForeSlicesWeighted(:, sc) = mean(fq.TC.ForeSlices(:, sc)) + ...
        min(1.5, nslc * sum(sum(fmrifg(:, :, sc))) / nforeg) .* ...
        (fq.TC.ForeSlices(:, sc) - mean(fq.TC.ForeSlices(:, sc)));
end
fq.TC.ForeSlicesWeighted(isnan(fq.TC.ForeSlicesWeighted)) = 0;

% estimate outliers of foreground (global)
[trashvar, outlest] = windsorize(fq.TC.Foreground, 5, 5);
fq.Quality.Outliers.Volumes = fq.Quality.Outliers.Volumes + double(outlest(:));
fq.Quality.Outliers.VolumesForeTC = find(outlest);
dvarin = abs(diff(fq.TC.Foreground(:)));
dvarin = [dvarin; 0] + [0; dvarin];
[trashvar, outlest] = windsorize(dvarin, 5, 5);
fq.Quality.Outliers.Volumes = fq.Quality.Outliers.Volumes + double(outlest(:));
fq.Quality.Outliers.VolumesForeTCDiff = find(outlest);

% estimate outliers of foreground (slices)
mdistest = madistd(fq.TC.ForeSlicesWeighted(:, var(fq.TC.ForeSlicesWeighted) > 0));
[trashvar, outlest] = windsorize(mdistest, 5, 5);
fq.Quality.Outliers.Volumes = fq.Quality.Outliers.Volumes + double(outlest(:));
fq.Quality.Outliers.VolumesForeSlicesTC = find(outlest);
dvarin = abs(diff(fq.TC.ForeSlicesWeighted));
dvarin = [dvarin; zeros(1, nslc)] + [zeros(1, nslc); dvarin];
mdistest = madistd(dvarin(:, var(dvarin) > 0));
[trashvar, outlest] = windsorize(mdistest, 5, 5);
fq.Quality.Outliers.Volumes = fq.Quality.Outliers.Volumes + double(outlest(:));
fq.Quality.Outliers.VolumesForeSlicesWDiff = find(outlest);

% get background noise sample
fmrins = sort(fq.Raw.StdevImage(fmribg));

% and compute mean over [0.25 ... 0.75] interval
fmrins = mean(fmrins(ceil(numel(fmrins) / 4):floor(3 * numel(fmrins) / 4)));

% compute signal to noise images
fq.Quality.GlobalSNRImage = (1 / fmrins) .* fq.Raw.MeanImage;
fq.Quality.LocalSNRImage = fq.Raw.MeanImage ./ fq.Raw.StdevImage;

% build temporal filtering argument
if opts.tempfset
    tfdct = opts.tempffrq;
    tfscf = 0;
else
    tfdct = Inf;
    tfscf = floor(nvol / opts.tempffrq);
end
tfstr = struct('tdim', 4, 'tempdct', tfdct, 'tempsc', tfscf, 'trobust', true);

% preliminary filter data
[fmrifilt, fx, fmrioutl] = tempfilter(fmridata, tfstr);
fq.TempFiltered.StdevImage = std(fmrifilt, [], 4);
fmrifns = sort(fq.TempFiltered.StdevImage(fmribg));
fmrifns = mean(fmrifns(ceil(numel(fmrifns) / 4):floor(3 * numel(fmrifns) / 4)));

% compute filtered signal to noise images
fq.Quality.TF_GlobalSNRImage = (1 / fmrifns) .* fq.Raw.MeanImage;
fq.Quality.TF_LocalSNRImage = fq.Raw.MeanImage ./ fq.TempFiltered.StdevImage;

% get filtered time courses
fq.TC.TF_Global = lsqueeze(mean(mean(mean(fmrifilt))));
fq.TC.TF_Slices = squeeze(mean(mean(fmrifilt)))';
fq.TC.TF_Foreground = fq.TC.TF_Global;
fq.TC.TF_ForeSlices = fq.TC.TF_Slices;
for vc = 1:nvol
    svol = fmrifilt(:, :, :, vc);
    fq.TC.TF_Foreground(vc) = mean(svol(fmrifg));
    for sc = 1:nslc
        sslc = svol(:, :, sc);
        fq.TC.TF_ForeSlices(vc, sc) = mean(sslc(fmrifg(:, :, sc)));
    end
end
fq.TC.TF_ForeSlices(isnan(fq.TC.TF_ForeSlices)) = 0;

% weigh slice timecourses by number of voxels
fq.TC.TF_ForeSlicesWeighted = fq.TC.TF_ForeSlices;
for sc = 1:nslc
    fq.TC.TF_ForeSlicesWeighted(:, sc) = mean(fq.TC.TF_ForeSlices(:, sc)) + ...
        min(1, nslc * sum(sum(fmrifg(:, :, sc))) / nforeg) .* ...
        (fq.TC.TF_ForeSlices(:, sc) - mean(fq.TC.TF_ForeSlices(:, sc)));
end
fq.TC.TF_ForeSlicesWeighted(isnan(fq.TC.TF_ForeSlicesWeighted)) = 0;

% estimate more outliers ...
[trashvar, outlest] = windsorize(fq.TC.TF_Foreground, 4, 10);
fq.Quality.Outliers.Volumes = fq.Quality.Outliers.Volumes + double(outlest(:));
fq.Quality.Outliers.VolumesTF_ForeTC = find(outlest);
mdistest = madistd(fq.TC.TF_ForeSlicesWeighted(:, var(fq.TC.TF_ForeSlicesWeighted) > 0));
[trashvar, outlest] = windsorize(mdistest, 4, 10);
fq.Quality.Outliers.Volumes = fq.Quality.Outliers.Volumes + double(outlest(:));
fq.Quality.Outliers.VolumesTF_ForeSlicesTC = find(outlest);

% get outlier time course
fq.TC.Outliers = fq.TC.Global;
fq.TC.OutlierSlices = fq.TC.Slices;
fgvox = 1 / sum(fmrifg(:));
fgsvox = ones(1, nslc);
for sc = 1:nslc
    fgsvox(sc) = 1 / sum(sum(fmrifg(:, :, sc)));
end
for vc = 1:nvol
    svol = fmrioutl(:, :, :, vc);
    fq.TC.Outliers(vc) = fgvox * sum(svol(fmrifg));
    for sc = 1:nslc
        sslc = svol(:, :, sc);
        fq.TC.OutlierSlices(vc, sc) = fgsvox(sc) * sum(sslc(fmrifg(:, :, sc)));
    end
end

% compute how much variance was explained by filtering
[fmrifc, fmrifr] = cov_nd(fmridata, fmrifilt);
fmrifr(fmrifr < 0) = 0;
fq.Quality.TP_SumSquares = 1 - sqrt(fmrifr);

% create corrected time series (for further computation)
fmrinflt = fmrifilt;

% set values with robust weight < 0.75 to NaN
fmrinflt(fmrioutl < 0.75) = NaN;

% and compute mean without NaNs
[fmrinflm, fmrinfge] = meannoinfnan(fmrinflt, 4);

% as well as number of "OK" volumes
fmrinfge = sum(fmrinfge, 4);

% replace NaNs with mean (without NaNs)
for vc = 1:nvol
    svol = fmrinflt(:, :, :, vc);
    svol(isnan(svol)) = fmrinflm(isnan(svol));
    fmrinflt(:, :, :, vc) = svol;
end

% now re-compute std, max and min of corrected time series
fmrinfst = sqrt((nvol - 1) .* (var(fmrinflt, [], 4) ./ (fmrinfge - 1)));
fq.Quality.TF_NoOutliers_StdevImage = fmrinfst;

% temporal signal (maximal contrast)
fmrinfmm = max(fmrinflt, [], 4) - min(fmrinflt, [], 4);

% give quality measure in numbers of BG noise
fq.Quality.TF_ContrastToNoise = (1 / fmrins) .* max(fmrinfmm - fmrinfst, 0);

% compute z-transformed filtered data time courses
fmrifilt = abs(ztrans(fmrifilt, 4));
fmrifilt(isnan(fmrifilt)) = 0;
fq.TC.TF_zScoreGlobal = fq.TC.Global;
fq.TC.TF_zScoreSlices = fq.TC.Slices;
for vc = 1:nvol
    svol = fmrifilt(:, :, :, vc);
    fq.TC.TF_zScoreGlobal(vc) = fgvox * sum(svol(fmrifg));
    for sc = 1:nslc
        sslc = svol(:, :, sc);
        fq.TC.TF_zScoreSlices(vc, sc) = fgsvox(sc) * sum(sslc(fmrifg(:, :, sc)));
    end
end

% compute outlier image
fmrioutl = (1 / nvol) .* (nvol - sum(fmrioutl .* fmrioutl, 4));
fq.Quality.Outliers.WeightImage = fmrioutl;

% some global quality time courses
fq.TC.Quality = ztrans([ ...
    fq.TC.Global, fq.TC.Foreground, fq.TC.Outliers, fq.TC.TF_zScoreGlobal, ...
    fq.TC.SmoothEst]);

% do not perform extended (moco) analyses?
if ~opts.motcor
    if opts.qasheet
        fmriqasheet(fq);
    end
    return;
end

% perform motion correction
trf = eye(4);
trf([1, 6, 11]) = opts.res;
trf(13:15) = 0.5 * (- (opts.res .* (1 + fmrisize(1:3)))');
[tfm, trp, fmridata] = ...
    rbalign(fmridata(:, :, :, 1), fmridata, struct('trfv1', trf, 'trfv2', trf));

% rescale motion parameters (for plottability)
fq.MotCorr.Params = zeros(nvol, 6);
for vc = 1:nvol
 	trpc = spmitrf(trp(:, :, vc));
    fq.MotCorr.Params(vc, :) = [trpc{1}, (180 / pi) * trpc{2}];
end

% create motion parameters nuisance regressors
mcnreg = [ ...
    ztrans(fq.MotCorr.Params), ...
    ztrans(fq.MotCorr.Params .* fq.MotCorr.Params)];

% compute some basic images
fq.MotCorr.MeanImage = (1 / nvol) .* sum(fmridata, 4);
fq.MotCorr.StdevImage = std(fmridata, [], 4);

% get voxels that are likely to be background
mfmrimean = mean(fq.MotCorr.MeanImage(:));
mdfmrimean = median(fq.MotCorr.MeanImage(fq.MotCorr.MeanImage < mfmrimean));
fmribg = (fq.MotCorr.MeanImage < mdfmrimean);
fmrifg = (fq.MotCorr.MeanImage > mfmrimean);
fmrifg = fmrifg & dilate3d(erode3d(fmrifg));

% and remove those with more than 5% (absolute) zero values (shifted out)
fmribg(sum(fmridata == 0, 4) > (0.05 * nvol)) = false;

% find biggest chunk
[cs, fmribg] = clustercoordsc(fmribg);
fmribg = (fmribg == maxpos(cs));
[cs, fmrifg] = clustercoordsc(fmrifg);
fmrifg = (fmrifg == maxpos(cs));
fq.MotCorr.Masks.Background = fmribg;
fq.MotCorr.Masks.Foreground = fmrifg;

% get some time courses
fq.MotCorr.TC.Global = lsqueeze(mean(mean(mean(fmridata))));
fq.MotCorr.TC.Slices = squeeze(mean(mean(fmridata)))';
fq.MotCorr.TC.Foreground = fq.MotCorr.TC.Global;
fq.MotCorr.TC.ForeSlices = fq.MotCorr.TC.Slices;
for vc = 1:nvol
    svol = fmridata(:, :, :, vc);
    fq.MotCorr.TC.Foreground(vc) = mean(svol(fmrifg));
    for sc = 1:nslc
        sslc = svol(:, :, sc);
        fq.MotCorr.TC.ForeSlices(vc, sc) = mean(sslc(fmrifg(:, :, sc)));
    end
end

% get average background noise sample
fmrins = sort(fq.MotCorr.StdevImage(fmribg));
fmrins = mean(fmrins(ceil(numel(fmrins) / 4):floor(3 * numel(fmrins) / 4)));

% compute signal to noise images
fq.MotCorr.Quality.GlobalSNRImage = (1 / fmrins) .* fq.MotCorr.MeanImage;
fq.MotCorr.Quality.LocalSNRImage = fq.MotCorr.MeanImage ./ fq.MotCorr.StdevImage;

% re-filter data (including motion correction parameters
fx = [fx, mcnreg, ones(size(fx, 1), 1)];
fmridata = reshape(fmridata, prod(fmrisize(1:3)), nvol);
[fmrifilt, fmrioutl] = fitrobustbisquare_img(fx, fmridata);
fmrifilt(:, end) = 0;
fmrifilt = fmridata - fmrifilt * fx';
fmridata = reshape(fmridata, fmrisize);
fmrifilt = reshape(fmrifilt, fmrisize);
fmrioutl = reshape(fmrioutl, fmrisize);

% get outlier time course
fq.MotCorr.TC.Outliers = fq.MotCorr.TC.Global;
fq.MotCorr.TC.OutlierSlices = fq.MotCorr.TC.Slices;
fgvox = 1 / sum(fmrifg(:));
fgsvox = ones(1, nslc);
for sc = 1:nslc
    fgsvox(sc) = 1 / sum(sum(fmrifg(:, :, sc)));
end
for vc = 1:nvol
    svol = fmrioutl(:, :, :, vc);
    fq.MotCorr.TC.Outliers(vc) = fgvox * sum(svol(fmrifg));
    for sc = 1:nslc
        sslc = svol(:, :, sc);
        fq.MotCorr.TC.OutlierSlices(vc, sc) = fgsvox(sc) * sum(sslc(fmrifg(:, :, sc)));
    end
end
[trashvar, outlest] = windsorize(fq.MotCorr.TC.Outliers, 3, 15);
fq.Quality.Outliers.Volumes = fq.Quality.Outliers.Volumes + double(outlest(:));
fq.Quality.Outliers.VolumesOutliersTC = find(outlest);
mdistest = madistd(fq.MotCorr.TC.OutlierSlices(:, var(fq.MotCorr.TC.OutlierSlices) > 0));
[trashvar, outlest] = windsorize(mdistest, 3, 15);
fq.Quality.Outliers.Volumes = fq.Quality.Outliers.Volumes + double(outlest(:));
fq.Quality.Outliers.VolumesOutlierSlicesTC = find(outlest);

% compute how much variance was explained by filtering
[fmrifc, fmrifr] = cov_nd(fmridata, fmrifilt);
fmrifr(fmrifr < 0) = 0;
fq.MotCorr.Quality.TP_SumSquares = 1 - sqrt(fmrifr);

% compute z-transformed filtered data time courses
fmrifilt = abs(ztrans(fmrifilt, 4));
fmrifilt(isnan(fmrifilt)) = 0;
fq.MotCorr.TC.TF_zScoreGlobal = fq.TC.Global;
fq.MotCorr.TC.TF_zScoreSlices = fq.TC.Slices;
for vc = 1:nvol
    svol = fmrifilt(:, :, :, vc);
    fq.MotCorr.TC.TF_zScoreGlobal(vc) = fgvox * sum(svol(fmrifg));
    for sc = 1:nslc
        sslc = svol(:, :, sc);
        fq.MotCorr.TC.TF_zScoreSlices(vc, sc) = fgsvox(sc) * sum(sslc(fmrifg(:, :, sc)));
    end
end
[trashvar, outlest] = windsorize(fq.MotCorr.TC.TF_zScoreGlobal, 3, 15);
fq.Quality.Outliers.Volumes = fq.Quality.Outliers.Volumes + double(outlest(:));
fq.Quality.Outliers.VolumesTF_Foreground_zTC = find(outlest);
mdistest = madistd(fq.MotCorr.TC.TF_zScoreSlices(:, var(fq.MotCorr.TC.TF_zScoreSlices) > 0));
[trashvar, outlest] = windsorize(mdistest, 3, 15);
fq.Quality.Outliers.Volumes = fq.Quality.Outliers.Volumes + double(outlest(:));
fq.Quality.Outliers.VolumeRatio = sum(fq.Quality.Outliers.Volumes > 2) / nvol;
fq.Quality.Outliers.VolumesTF_ForeSlices_zTC = find(outlest);

% compute outlier image
fmrioutl = (1 / nvol) .* (nvol - sum(fmrioutl .* fmrioutl, 4));
fq.Quality.Outliers.MC_WeightImage = fmrioutl;

% output sheet ?
if opts.qasheet
    fmriqasheet(fq);
end
